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Characterization of small RNAs in microsporidian 


Nosema bombycis 
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Abstract; [Aim] The microsporidian Nosema bombycis can produce epizootic outbreaks in populations of the 
domestic silkworm (Bombyx mori) by the major mode of vertical transmission and horizontal transmission. In 
order to explore the genomic defense system against repetitive elements and potential means of gene regulation in 
N. bombycis, we conducted the survey of small RNAs associated with transposons and identification of the 
potential miRNAs in this silkworm parasite. { Methods] Total RNAs were extracted from silkworm midguts 
infected by N. bombycis spores and small RNA cDNA libraries were constructed for deep sequencing. After 
Solexa sequencing, we performed the functional characterization of N. bombycis small RNAs based on the 
reference genomic sequence through bioinformatics. [ Results] N. bombycis has predominantly small RNAs with 
the length of 24 and 25 nt, most of which show a strong base bias for U at the 5’ end. Redundant small RNAs 
were associated with transposons in N. bombycis, and the amount of antisense small RNAs aligned with 
transposons was much more than that of sense small RNAs. Thirty-one candidate miRNAs were characterized by 
means of computational prediction, and certain miRNAs were shared by other Nosema species, suggesting that 
they may evolve conservatively. [ Conclusion] There should exist small RNAs genomic defense system against 
transposons in N. bombycis. Tentative identification of potential thirty-one miRNAs provides a basis for further 
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functional research of their targets. 
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1 INTRODUCTION 


Since small RNAs play important roles in the 
development of individual organisms and regulation 
( Kawaoka et al., 2008). In recent decades, 
numerous studies have focused on small non-coding 
RNAs from various organisms. According to the 
mechanism of their biosynthesis pathways, small 
RNAs can be divided into three different categories : 
microRNAs (miRNAs), small interfering RNAs 
(siRNAs) and PIWI-interacting RNAs ( piRNAs) 
(Lau et al., 2009). siRNAs were generated as 21 - 
25 nt in length through a long double-stranded RNA 
molecule digested by a Dicer enzyme, and then 
silenced the transcription of targeted genes by 
binding with AGO proteins (Ender et al., 2008; 
Hick and Meister, 2008; Obbard and Finnegan, 
2008 ). Mature miRNAs, about 19 - 25 nt in 
length, were generated from the precursor miRNAs 
(premary miRNAs ) containing hairpin structure 
through the action of the RNase III family enzymes 














Drosha and Dicer, and were then bound by Ago- 
subfamily proteins ( Sinkkonen et al., 2008; Kim et 
al., 2009). miRNA regulated gene expression at the 
post-transcriptional level, and mediated several 
important biochemical pathways ( MacRae et al., 
2007). piRNAs, initially known as rasiRNAs, were 
a class of newly discovered small RNAs with the 
length of 26 -32 nt (Girard et al., 2006; Houwing 
et al., 2007; Wang et al., 2014). Generation of 
piRNAs depended on the role of PIWI subfamily 
proteins, but not Dicer enzyme ( Vagin et al., 
2006). In Drosophila, piRNAs were initially named 
RNAs 


(rasiRNAs ), due to mainly derivation from the 


as ` repeat-associated small ` interfering 
repetitive DNA sequences ( Saito et al., 2006; 
Brennecke et al., 2007). 

The function of small RNAs had been reported 
in several fungi and protists. For instance, small 
RNAs of Candida albicans were able to silence 
endogenous retrotransposon activity ( Drinnenberg et 


al., 2009); small RNAs from Schizosaccharomyces 

















tt 








Lamp, AC DOR OTF ACR RM H (2013 AA102507) ; DS 














SERGE SH Al (31001037 31270138, 31402142) ; ARAL BI Bl 























AYIR (201410637005 ) ; EET AMES RUS WEIN ILRI H (cste2015 jcyjA80020 , cste2015 jcyjA80014 ) 
EET: En w, 1991 469 AE, BLA, wt. Mä the: 9 BE , E-mail: pql19910907@ 126. com 








* HTHVES Corresponding author, E-mail; xujinshan2003@ aliyun. com 
WAR A H Received; 2015-06-16; 25% H HH Accepted; 2015-09-29 








1214 EE} Acta Entomologica Sinica 58 ZS 





pombe can regulate heterochromatin assembly and 
maintain heterochromatin silencing ( Bühler and 
Moazed, 2007 ). 


candidate miRNAs have been screened and identified 


In Entamoeba histolytica, 17 


by using computational methods (De et al., 2006). 
Microsporidia were obligated intracellular parasites 
invading vertebrates and invertebrates. Molecular 
evidence proved that microsporidia were closely 
related to fungi (Lee et al., 2008). Up to date, at 
least twenty-one genomes of microsporidia species 
have been sequenced ( Parisot et al., 2014; 
Desjardins et al., 2015; Pombert et al., 2015). 
However, no thorough survey of small RNAs in any 
of these species mentioned above was reported. 
Nosema bombycis , a typical species of microsporidia , 
can infect the silkworm ( Bombyx mori) through 
vertical transmission from the mother host to their 
progenitive eggs and horizontal transmission via 
mouth, and consequently, has a devastating impact 
on the silkworm industry (Ma et al., 2013; Zhao et 
al., 2015). The genome of N. bombycis is 15. 3 
Mbp in size, 40% 
repetitive elements (Pan et al., 2013). However, 


of which composes of the 


the defense mechanism against the amplification of 
repetitive elements remains unclear in the N. 
bombycis genome. It had been reported that small 
RNAs could defense against transposons in the 
silkworm ( Kawaoka et al., 2008 ). 
whether small RNAs are associated with transposon 


Therefore , 


defense in N. bombycis would be worth excavating. 
In the study, we conducted the survey of N. 
bombycis small RNAs to identify the small RNAs 
associated with transposons and the potential 
miRNAs in N. bombycis, which provides a basis for 
further study of the biological function of N. 


bombycis small RNAs. 


2 MATERIALS AND METHODS 


2.1 Test insect 

B. mori strain Dazao reserved in our laboratory 
was cultivated as the 3rd instar larvae after eggs 
hatched, and then the 3rd instar larvae were infected 
by spores of N. bombycis CQ1, which were isolated 
China, 
conserved in China Veterinary Culture Collection 
Center (CVCC number: 102059). At 4 d post 


infection, some tissues were collected to assay the 


from infected silkworms in Chongqing, 


condition of spores propagation through observation 
under microscope. Finally, whole bodies of day-3 
Sth instar larvae were grinded and the total number 
of spores was counted by using hemocytometer. 
2.2 Genomic data resources 

Genomic sequences and gene annotation of N. 


bombycis CQ1 were downloaded from public Silkworm 
Pathogen Database (SilkPathDB, http://silkpathdb. 
swu. edu. cn/silkpathdb/ftpserver ) and GenBank 
(accession no. ACJZO1000001 — ACJZ01003558 ) . 
2.3 Construction of small RNA cDNA libraries 
for deep sequencing 

Whole bodies of day-3 5th instar larvae of B. 


mori infected with total 10° spores of N. bombycis 


CQ1 were collected for RNA isolation. After 
purification using Trizol, total RNAs were 
immediately stored in ethanol and at ` — 80° for 


further use. The total RNA was separated by 15% 
polyacrylamide gel electrophoresis, and recovered 
about 16 - 30 nt small RNAs. The 5’ RNA adapter 
( 5'-GUUCAGAGUUCUACAGUCCGACGAUC-3’) was 
ligated to the RNAs pool with T4 RNA ligase. The 3’ 
RNA adapter (5’-pUCGUAUGCCGUCUUCUGCUU 
GidT-3') was subsequently ligated to precipitated 
RNAs using T4 RNA ligase. Small RNAs ligated 
with adaptors were then subjected to RT-PCR to 
produce sequencing libraries. PCR products were 
purified and small RNA libraries were sequenced 
using Solexa, a massively parallel sequencing 
technology. After removing adaptor sequences and 
low quality (Q < 30) sequences, the clean reads 
were reserved. 
2.4 Classification and functional annotation of 
N. bombycis small RNAs 

Raw clean small RNA reads were aligned with 
sequences by Bowtie 


(mismatch <2 bp) software. Only those matched to 


genomic 


N. bombycis 


genomic sequence were retrieved and called as N. 
bombycis small RNAs. If small RNAs gave hits for 
repetitive sequence of N. bombycis at least once, 
these sRNAs were categorized as repeat-associated 
small interfering RNAs ( rasiRNAs). 


gave hits for non-coding sequence region, they were 


If sequences 


classified as ncRNAs. Likewise, those that gave hits 
for protein-coding regions were classified as ‘ gene’ . 
Others that matched the rDNAs or tRNAs sequence 
were classified as ‘rDNAs’ or ‘tRNAs’ , respectively. 
Those ncRNAs of N. bombycis were picked up to 
predict potential miRNAs by MIREAP ( https:// 
sourceforge. net/projects/mireap/). Secondary RNA 
structures of candidate miRNAs were visualized using 
RNAfold ( Hofacker and Stadler, 2006). Three 
miRNA databases, including MiRAlign ( http:// 
bioinfo. au. tsinghua. edu. cn/miralign/miralign _ 
dscems. htm), MiRbase (http://microrna. sanger. 
and Rfam ( http;//rfam. 


sanger. ac. uk/), were retrieved for homologous 


ac. uk/sequences/ ) 


search of potential miRNAs identified in this study. 
To intuitively observe homologous sequences of N. 
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miRNAs in other related 
both mature miRNAs and 
precursor miRNAs from N. bombycis were Blast- 


bombycis candidate 


microsporidia genomes, 


aligned with the public genomic sequences of six 
microsporidia species including Nosema antheraeae , 
Nosema cuniculi , 


ceranae , Encephalitozoon 


Encephalitozoon intestinalis, Enterocytozoon bieneusi 
and Antonospora locustae. The parameters for the 
BLAST search were set as following; matching region 
covered more than 60% of whole query sequence and 
pairwise identity thresholds were 80% at nucleic 
acid level. Subsequently, statistical analysis was 


performed to investigate the evolutionary conservation 
of these candidate miRNAs. 


3 RESULTS 


3.1 Sequencing and classification of small 
RNAs from N. bombycis 

A total of 3 844 786 raw reads were obtained by 
Solexa sequencing. After trimming low-quality, poly 
(N) adapter and adapter sequences, 2 367 321 
clean reads of small RNAs remained. Through 
matching to N. bombycis genomic sequences, 661 609 
reads were collected as small RNAs sequences of N. 
bombycis ( Table 1). It was notable that there were 
still 1 705 712 unclassified reads expect N. bombycis 
and these reads should be 
originated from host domestic silkworm ( B. mori), 
because the mixed RNAs isolated from midguts of the 


small RNAs sequences, 
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domestic silkworm infected by N. bombycis must 
contain massive RNAs of the silkworm in this study. 
N. bombycis small RNAs showed a unimodal length 
distribution, with a peak at 25 nt (Fig. 1: A). 
Seventy percent of cloned N. bombycis small RNAs 
contained a 5’ uridine residue (Fig. 1: B), 
indicative of bases bias in the 5’ terminus. 

Small RNAs sequences obtained were classified 
by matching N. bombycis genomic annotated 
database (Table 1). Of the total 661 609 small 
RNAs, 460 916 reads (75. 93% ) matched to the 
non-coding regions of the genome (ncRNA) , 162 691 
reads (21.68% ) matched to the repetitive sequences 
of the genome (rasiRNA), and 1. 73%, 0. 65% 
and 0.01% of small RNAs matched to the protein- 
coding regions, tRNAs and rDNAs, respectively. 
All these indicated that ncRNAs and 


rasiRNAs were mostly predominant in all N. 


results 


bombycis small RNAs sequences. 


Table 1 Sequencing data of Nosema bombycis small RNAs 


Total number of 
Number of reads ; ` 
nucleic acids 








Raw reads 3 844 786 13 567 510 
Low quality ( < Q30) 1 418 249 49 638 715 
Poly( N) reads 7 870 275 450 
Clean reads ( adapter trimming) 2 367 321 61 598 363 
Matching N. bombycis 661 609 15 997 478 
Unmatched 1705712 45600885 ` 1 705 712 45 600 885 
B 100 
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Fig. 1 


Characterization of Nosema bombycis small RNAs 


A; Length distributions of sequenced small RNA reads; B; Nucleotide compositions of sequenced small RNAs; C: Proportions of annotated small 


RNAs in N. bombycis. 5'1 —5; 1st —5th bases at the 5’ terminus; 3'1 -5: 


Ist — 5th bases at the 3’ terminus. 
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3.2 Small RNAs matched with transposons of 
N. bombycis 

Most of N. bombycis rasiRNAs were found to be 
annotated N. bombycis 
transposable elements, including LTR, non-LTR and 
DNA types. rasiRNAs corresponding to various types 
of annotated transposons were summarized in Table 
2. Of the total 162 691 rasiRNAs, 8 417 rasiRNAs 
matched with LTR retrotransposons and 6 516 
rasiRNAs matched with DNA transposons. 
Interestingly, rasiRNAs gave hits for all the 
annotated V. bombycis LTR retrotransposons expect 
for the Nbr3 and Nbr8 families. Generally, N. 
bombycis rasiRNAs gave more frequent hits toward 
antisense orientation of transposons than those hits 
toward sense orientation (Table 2). Two full-length 
LTR retrotransposons ( Nbr4 and Nbr7), as well as 


associated with well 


one DNA transposon (NbTcl ) , were then chosen to 
analyze the distribution of sRNAs in the complete 
elements. After plotting all of the rasiRNAs that 
were mapped to the targets, the characteristic 
phasing patterns were identified as in Fig. 2. 
rasiRNAs are distributed inconsecutively in the three 
complete elements and give much more frequent hits 
for antisense orientation than for sense orientation. 
In the case of Nbr4, rasiRNA is distributed 
throughout the element, and two peaks of antisense 
rasiRNAs were found at the sites of 1 200 and 4 100 
nt, respectively (Fig. 2: A). In the case of Nbr7, 
rasiRNA mapped all the regions excluding the two 
termini, and a peak of antisense rasiRNAs was 
present at site of 1 300 nt (Fig. 2: B). For NbTcl, 
sense and antisense rasiRNAs are mainly distributed 
in the middle of the element (Fig. 2; C). 


Table 2 Number and percent of small RNAs corresponding to the different types of annotated 
transposons in Nosema bombycis 














Type Family Sense Antisense Total Percent References 
MITE 1 561 1 012 2513 14.8 He et al., 2015 
DNA Tcl /mariner 750 2 194 2 944 16.9 Xu et al., 2010 
Tiger-like 305 694 999 5.7 Pan et al., 2013 
Nbr1 -Nbr8 869 1 663 2 532 14.5 Xu et al., 2006 
LTR Nbr9-Nbr14 636 1 520 2 156 12.4 Xiang et al., 2010 
Others 1 942 1 787 3 729 21.4 Pan et al., 2013 
Non-LTR NbR4 576 1914 2 490 14.3 Liu et al., 2014 
Total 6 639 10 784 17 423 
3.3 Identification of N. bombycis potential miRNAs m0003-5p and NB-m0004-5p, were clustered 


Since 76% of small RNA sequences were 
ncRNA sequences, we intended to explore whether 
there were potential miRNAs derived from ncRNAs 
in N. bombycis. After all the neRNAs were mapped 
to genome loci, only those 300 bp genomic regions 
harboring more than 50 ncRNA tags, were chosen as 
the potential sites to produce miRNA candidate. 
More than 1 000 sites were picked up to predict 
miRNAs through MIREAP software under the 
criterion of energy < — 20 kal/mol. And then 
“precursor miRNAs’ containing intricate stem-loop 
structures were discarded manually by means of 
RNAfold ( Hofacker and Stadler, 2006 ). 
Ultimately, 31 candidate miRNAs and miRNA 
precursors were generated from the N. bombycis 
genome loci sequences. The characteristics of these 
31 potential miRNAs were shown in Table 3, and 
several cases of miRNA stem-loop structure were 
visualized in Fig. 3( A). While 3 kb interval was 
defined as the threshold value between two miRNAs 
( Altuvia et al., 2005), several candidate miRNA 
genes identified here were arrayed in cluster in the 
N. bombycis genome. Two miRNAs, named as NB- 


together since there was only 530 bp interval between 
them (Fig. 3: B). In order to investigate the 
presence and conservation of potential miRNAs in 
other microsporidium species, 31 mature miRNAs 
and precursor miRNAs from N. bombycis were Blast- 
aligned with the public genomic data of six 
microsporidia species, which were N. antheraeae, 
N. ceranae, E. cuniculi, E. intestinalis, E. 
bieneusi and A. locustae, respectively. As shown in 
Fig. 3 (C), it was found that not only 14 mature 
miRNAs but also five precursor miRNAs from N. 
bombycis share sequence similarity in N. antheraeae. 
However, only three homologs of N. bombycis mature 
miRNAs exist in N. ceranae or E. bieneusi, and no 
hit was found in E. cuniculi, E. intestinalis and A. 
locustae. These results implied that miRNAs might 
exist in the common ancestor of microsporidia and 
undergo the loss events during microsporidia genomic 
of miRNA targets 
conducted tentatively through alignment of miRNAs 
with the 3’ untranslated region (UTR) of any N. 
bombycis genes. As a result, one candidate miRNA 


named as NB-m0189-3p matched perfectly with the 


evolution. Prediction were 
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Fig. 2 Density plot of rasiRNAs assigned to Nbr3 (A), Nbr7 (B) and NbTcl (C) 


Red curve; Sense; Blue curve; Antisense. 


3'-UTR of NBO-1020036 gene (Fig. 3: D). 
4 DISCUSSION 
Although the key genes responsible for RNAi 


mechanism have been found recently in microsporidia 
species (Paldi et al., 2010; Heinz et al., 2012; Wang 
et al., 2015), there was no report about the small 
RNAs in microsporidia. In this study, we performed a 
systematic and integrative analysis of small RNAs in 
the microsporidia typical species, silkworm parasite N. 
bombycis. After Solexa sequencing and classification of 
N. bombycis small RNAs, 75.93% and 21.68% were 
found to be ncRNAs and rasiRNAs, respectively. 
Notably, strong U bias at 5’ ends of N. bombycis small 
RNAs was displayed. This phenomenon was therein 


consistent with those small RNAs reported in the 
domestic silkworm ( Kawaoka et al., 2008). Our 
results showed that lots of rasiRNAs gave hits to most 
well-annotated transposons and mapped their targets in 
differential manners, indicating that they should have 
the role in regulating activity of transposon. Moreover, 
the strong strand bias of N. bombycis rasiRNAs toward 
antisense orientation was usually observed in our study , 
which might be due to nucleotide constitution of both 
sense and antisense strands. Previous studies reported 
that more than 20% of N. bombycis genome sequences 
were comprised of transposable elements, which have 
an important impact on genomic and gene organization 
(Xu et al., 2010; Pan et al., 2013), and several 


transposons still have transcribed activity in N. bombycis 
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Fig. 3 Characterization of Nosema bombycis candidate miRNAs 
A; Visualization of secondary structure of four candidate miRNAs from N. bombycis; B; miRNA clusters in the genome; C; Homologous sequences of 


N. bombycis miRNAs present in different species; D; Alignment of miRNAs to the 3’ UTRs of potential targets. 


(Xiang et al., 2010; Xu et al., 2010). Therefore, so 
in N. bombycis 
genome, it raised the question if there is molecular 


many transposons were present 
mechanism of inhibiting the mobilization of transposon 
and therein keeping the genome stability. Based on our 
bioinformatic finding in this study, it can imply that 
these rasiRNAs should be involved in genomic defense 
system against transposons in N. bombycis. However, 
the main question remains to be unresolved was that 
how these rasiRNA can silence the transposon in N. 
bombycis genome. Some studies have reported that 
piRNAs cluster produced from a Flamenco locus in X 
chromosome could mediate the transcriptional level of 
three retrotransposons gypsy, ZAM and Idefix in 
Drosophila melanogaster ( Brennecke et al., 2007) , 
and an interaction between Aubergine-associated 
antisense piRNAs and AGO3-associated sense piRNAs 
was regarded to lead to the efficient silencing of 
transposons in Drosophila (Gunawardane et al., 2007). 
Those findings mentioned above offered the clue to see 
whether putative AGO gene can interact with rasiRNA 
and then trigger the silence of active transposon of N. 
bombycis. 

In this study, structural 
features, including hairpin structure, the length of 


we found several 
stem-loop, bulge size and thermodynamic stability, 
which were usually considered as the standards for 


bioinformatic prediction of miRNAs (De et al., 2006; 
Chen et al., 2009). On this basis, 31 potential 


candidate miRNAs were predicted in N. bombycis, 
some of which were clustered in the certain genomic 
regions and then might be transcribed as polycistronic 
structure. In the amphioxus, 45 miRNAs constituted 
17 compact clusters, some of which were even 
conserved in vertebrates (Chen et al., 2009). So our 
results indicated that the existence of miRNAs clusters 
might be in most organisms harboring 
miRNAs. However, the roles of these candidate 
miRNAs of N. bombycis need to be further resolved. 
Even if the targeted gene NBO-10g0036 has been 
predicted by using NB-m0189-3p as query to match 
perfectly with the 3’-UTR of gene, more evidence is 
needed to determine whether the NB-m0189-3p can 
downregulate the transcriptional level of NBO-1020036 
gene based on experiment. Candidate miRNAs share 
sequence similarity among N. bombycis, N. antheraeae 
and N. ceranae, indicating that they should be present 


commonly in Nosema. No homolog of miRNAs exists in 


common 


E. cuniculi and E. intestinalis, which was consistent 
with their traits of genomic compaction. However, it 
was inexplicable that three homologs of miRNA were 
also found in the E. bieneusi, whose genome was 
considered to be extremely reduced ( Akiyoshi et al., 
2009). So, more evidence is needed to explain this 
observation. 

All together, this was the first time to identify 
small RNAs in N. bombycis, which should undoubtedly 


contribute to characterize small RNAs from other 
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In addition, the 
bombycis rasiRNAs and potential 


microsporidian species in future. 
findings of N. 
miRNAs will help greatly to clarify the origin and 
function of microsporidia small RNAs. 
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